/*
 * Project CHEN Rendering Engine : /geometry.hpp
 * Created: 2020-12-16, Last modified:  2021-02-16
 * 
 * This is the default descriptor.
 * 
 * Copyright (C) 2021 Hineven (hineven@pku.edu.cn)
 * 
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#ifndef _H_GEOMETRY
#define _H_GEOMETRY
#include "efloat.hpp"
#include "heads.hpp"
namespace chen {

inline bool isNan (const Vector2 & v) ;
inline bool isNan (const Vector3 & v) ;
inline bool isNan (const Point2 & v) ;
inline bool isNan (const Point3 & v) ;
inline bool isNan (const Normal & v) ;

class Vector2 {
	public:
    Float x, y;
    inline Vector2 (): x(0), y(0) {}
    inline explicit Vector2 (const Float & f): x(f),  y(f) {}
    inline Vector2 (Float a, Float b): x(a), y(b) {assert(!isNan(x) && !isNan(y));}
    inline explicit Vector2 (const Point2 & p) {*this = *((Vector2*)((void*)(&p)));}
    inline Vector2 operator + (const Vector2 & b) const {return Vector2(x+b.x, y+b.y);}
    inline Vector2 operator - (const Vector2 & b) const {return Vector2(x-b.x, y-b.y);}
    inline Vector2 operator * (const Float & b) const {assert(!isNan(b)); return Vector2(x*b, y*b);}
    inline Vector2 operator / (const Float & b) const {assert(b != 0); return Vector2(x/b, y/b);}
    inline Vector2 & operator += (const Vector2 & b) {x += b.x, y += b.y; return *this;}
    inline Vector2 & operator -= (const Vector2 & b) {x -= b.x, y -= b.y; return *this;}
    inline Vector2 & operator *= (const Float & b) {assert(!isNan(b)); x *= b, y *= b; return * this;}
    inline Vector2 & operator /= (const Float & b) {assert(b != 0); x /= b, y /= b; return * this;}
    inline Vector2 operator - () const {return Vector2(-x, -y);}
    inline Float lengthSquared () const {return x*x+y*y;}
    inline Vector2 normalized () const {
        Float fac = 1/std::sqrt(lengthSquared());
        assert(!isNan(fac));
        return (*this)*fac;
    }
    inline void resize (const Float & f) {
        Float fac = (f/std::sqrt(lengthSquared()));
        assert(!isNan(fac));
        *this = *this*fac;
    }
    inline Vector2 resized (const Float & f) const {
        Float fac = (f/std::sqrt(lengthSquared()));
        assert(!isNan(fac));
        return (*this)*fac;
    }
};
class Vector3 {
	public:
    Float x, y, z;
    inline Vector3 (): x(0), y(0), z(0) {}
    inline explicit Vector3 (const Float & f): x(f), y(f), z(f) {}
    inline Vector3 (Float a, Float b, Float c): x(a), y(b), z(c) {assert(!isNan(x) && !isNan(y) && !isNan(z));}
    inline explicit Vector3 (const Point3 & p) {*this = *((Vector3*)((void*)(&p)));}
    inline explicit operator Point3 () const ;
    inline explicit operator Normal () const ;
    inline Vector3 operator + (const Vector3 & b) const {return Vector3(x+b.x, y+b.y, z+b.z);}
    inline Vector3 operator - (const Vector3 & b) const {return Vector3(x-b.x, y-b.y, z-b.z);}
    inline Vector3 operator * (const Float & b) const {assert(!isNan(b)); return Vector3(x*b, y*b, z*b);}
    inline Vector3 operator / (const Float & b) const {assert(b != 0); return Vector3(x/b, y/b, z/b);}
    inline Vector3 & operator += (const Vector3 & b) {x += b.x, y += b.y, z += b.z; return *this;}
    inline Vector3 & operator -= (const Vector3 & b) {x -= b.x, y -= b.y, z -= b.z; return *this;}
    inline Vector3 & operator *= (const Float & b) {assert(!isNan(b)); x *= b, y *= b, z *= b; return * this;}
    inline Vector3 & operator /= (const Float & b) {assert(b != 0); x /= b, y /= b, z /= b; return * this;}
    inline Vector3 operator - () const {return Vector3(-x, -y, -z);}
    inline Float lengthSquared () const {return x*x+y*y+z*z;}
    inline Vector3 normalized () const {
        Float fac = 1/std::sqrt(lengthSquared());
        assert(!isNan(fac));
        return (*this)*fac;
    }
    inline Vector3 normalizedSafe () const {
        Float len = std::sqrt(lengthSquared());
        if(len == 0) return {0, 0, 0};
        Float fac = 1/len;
        return (*this)*fac;
    }
    inline void resize (const Float & f) {
        Float fac = (f/std::sqrt(lengthSquared()));
        assert(!isNan(fac));
        *this = (*this)*fac;
    }
    inline Vector3 resized (const Float & f) const {
        Float fac = (f/std::sqrt(lengthSquared()));
        assert(!isNan(fac));
        return (*this)*fac;
    }
    inline Float & operator [] (int idx) {
        if(idx == 0) return x;
        else return idx == 1 ? y : z;
    }
    inline const Float & operator [] (int idx) const {
        if(idx == 0) return x;
        else return idx == 1 ? y : z;
    }
    inline void permute (int a, int b, int c) {
        Float nx = a == 0 ? x : (a == 1 ? y : z);
        Float ny = b == 0 ? x : (b == 1 ? y : z);
        Float nz = c == 0 ? x : (c == 1 ? y : z);
        x = nx, y = ny, z = nz;
    }
};

inline Vector3 operator * (const Float & a, const Vector3 & b) {
    return Vector3(a*b.x, a*b.y, a*b.z);
}

template<class T> inline T operator + (const Vector3 & a, const T & b) {
    return Vector3(a.x+b.x, a.y+b.y, a.z+b.z);
}
template<class T> inline Float dot (const Vector3 & a, const T & b) {
    return a.x*b.x+a.y*b.y+a.z*b.z;
}
inline Vector2 abs (const Vector2 & a) {return Vector2(abs(a.x), abs(a.y));}
inline Vector3 abs (const Vector3 & a) {return Vector3(abs(a.x), abs(a.y), abs(a.z));}
inline Vector3 cross (const Vector3 & a, const Vector3 & b) {
    return Vector3(a.y*b.z-a.z*b.y, a.z*b.x-a.x*b.z, a.x*b.y-a.y*b.x);
}

// We assume that all normals are of unit length, except normal differentials
class Normal {
private:
    inline Float lengthSquared () const {return std::sqrt(x*x+y*y+z*z);}
public:
    Float x, y, z;
    inline Normal (): x(0), y(0), z(1) {}
    inline Normal (Float a, Float b, Float c): x(a), y(b), z(c) {}
    inline explicit Normal (const Point3 & n) {*this = *((Normal*)((void*)&n));}
    inline explicit Normal (const Vector3 & p) {*this = *((Normal*)((void*)&p));}
    inline explicit operator Vector3 () const {return *((Vector3*)((void*)this));}
    inline explicit operator Point3 () const ;
    inline Normal operator - () const {return Normal(-x, -y, -z);}
    inline Vector3 operator * (const Float & b) const {
        return Vector3(x*b, y*b, z*b);
    }
    inline Vector3 getPerpendicularUnit () const {
        if(!eequal(abs(x), 1)) return cross((Vector3)*this, Vector3(1, 0, 0)).normalized();
        else return cross((Vector3)*this, Vector3(0, 1, 0)).normalized();
    }
    inline Vector3 coord(const Normal & y, const Vector3 & v) const
    {
        return Vector3(dot(cross((Vector3)y, (Vector3)*this), v),  dot((Vector3)y, v), dot((Vector3)*this, v));
    }
    inline Vector3 rcoord(const Normal & y, const Vector3 & v) const
    {
        return y*v.y+(*this)*v.z+cross((Vector3)y, (Vector3)*this)*v.x;
    }
    inline Vector3 resized (Float l) const {
        return Vector3(x*l, y*l, z*l);
    }
    inline bool sideWith (const Vector3 & v) const {
        return x*v.x+y*v.y+z*v.z > (Float)0;
    }
    inline Float weigh (const Vector3 & v) const {
        return x*v.x+y*v.y+z*v.z;
    }
    inline void permute (int a, int b, int c) {
        Float nx = a == 0 ? x : (a == 1 ? y : z);
        Float ny = b == 0 ? x : (b == 1 ? y : z);
        Float nz = c == 0 ? x : (c == 1 ? y : z);
        x = nx, y = ny, z = nz;
    }
};

template<class T> inline Float dot (const Normal & a, const T & b) {
    return a.x*b.x+a.y*b.y+a.z*b.z;
}
class Point2 {
	public:
    Float x, y;
    Point2 ();
    explicit Point2 (const Float & f);
    Point2 (Float a, Float b);
    Point2 (const Point2 & v);
    explicit Point2 (const Vector2 & p) ;
    explicit operator Vector2 () const ;
    Point2 & operator = (const Point2 & b) ;
    Point2 operator + (const Vector2 & b) const ;
    Point2 operator - (const Vector2 & b) const ;
    Point2 operator * (const Float & b) const ;
    Point2 operator / (const Float & b) const ;
    Point2 operator *= (const Float & b) ;
    Point2 operator - () const ;
    Float lengthSquared () const ;
    Point2 normalized () const ;
    void resize (const Float & f) ;
    Point2 resized (const Float & f) const ;
};
class Point3 {
public:
    Float x, y, z;
    inline Point3 ():x(0), y(0), z(0) {}
    inline explicit Point3 (const Float & f): x(f), y(f), z(f) {}
    inline Point3 (Float a, Float b, Float c): x(a), y(b), z(c) {}
    inline explicit Point3 (const Normal & n): x(n.x), y(n.y), z(n.z) {}
    inline explicit Point3 (const Vector3 & p): x(p.x), y(p.y), z(p.z) {}
    inline explicit operator Normal () const {return *((Normal*)((void*)this));}
    inline explicit operator Vector3 () const {return *((Vector3*)((void*)this));}
    inline const Float & operator [] (int idx) const {
        if(idx == 0) return x;
        return idx == 1 ? y : z;
    }
    inline Float & operator [] (int idx)  {
        if(idx == 0) return x;
        return idx == 1 ? y : z;
    }
    inline Point3 operator + (const Vector3 & b) const {
        return Point3(x+b.x, y+b.y, z+b.z);
    }
    inline Point3 operator - (const Vector3 & b) const {
        return Point3(x-b.x, y-b.y, z-b.z);
    }
    inline Vector3 operator - (const Point3 & b) const {
        return Vector3(x-b.x, y-b.y, z-b.z);
    }
    inline Point3 operator * (const Float & b) const {
        return Point3(x*b, y*b, z*b);
    }
    inline Point3 operator - () const {
        return Point3(-x, -y, -z);
    }
    inline Point3 & operator += (const Vector3 & b) {
        x += b.x, y += b.y, z += b.z;
        return *this;
    }
    inline Point3 & operator -= (const Vector3 & b) {
        x -= b.x, y -= b.y, z -= b.z;
        return *this;
    }
    inline Point3 & operator *= (const Float & b) {
        assert(!isNan(b));
        x *= b, y *= b, z *= b;
        return *this;
    }
    inline Point3 & operator /= (const Float & b) {
        assert(!isNan(b) && b);
        x /= b, y /= b, z /= b;
        return *this;
    }
    inline Float lengthSquared () const {
        return x*x+y*y+z*z;
    }
    inline Point3 resized (Float l) const {
        Float f = l/std::sqrt(lengthSquared());
        assert(!isNan(f));
        return Point3(x*f, y*f, z*f);
    }
    inline Point3 normalized () const {
        Float f = 1/std::sqrt(lengthSquared());
        assert(!isNan(f));
        return Point3(x*f, y*f, z*f);
    }
    inline void permute (int a, int b, int c) {
        Float nx = a == 0 ? x : (a == 1 ? y : z);
        Float ny = b == 0 ? x : (b == 1 ? y : z);
        Float nz = c == 0 ? x : (c == 1 ? y : z);
        x = nx, y = ny, z = nz;
    }
};
template<class T> inline T operator + (const Point3 & a, const T & b) {
    return Point3(a.x+b.x, a.y+b.y, a.z+b.z);
}

bool isNan (const Vector2 & v) {return isNan(v.x) || isNan(v.y);}
bool isNan (const Vector3 & v) {return isNan(v.x) || isNan(v.y) || isNan(v.z);}
bool isNan (const Point2 & v) {return isNan(v.x) || isNan(v.y);}
bool isNan (const Point3 & v) {return isNan(v.x) || isNan(v.y) || isNan(v.z);}
bool isNan (const Normal & v) {return isNan(v.x) || isNan(v.y) || isNan(v.z);}

Vector3::operator Point3() const {return *((Point3*)((void*)this));}
Vector3::operator Normal() const {return *((Normal*)((void*)this));}
Normal::operator Point3() const {return *((Point3*)((void*)this));}

inline Point2 abs (const Point2 & a) {return Point2(abs(a.x), abs(a.y));}
inline Point3 abs (const Point3 & a) {return Point3(abs(a.x), abs(a.y), abs(a.z));}




Point3 offsetOrigin (const Point3 & o, const Vector3 & err, const Normal & n) ;

template<class T> inline Float dot (const Point3 & a, const T & b) {
    return a.x*b.x+a.y*b.y+a.z*b.z;
}

inline void coordinateSystem (const Normal & x, Vector3 * y, Vector3 * z) {
    *y = x.getPerpendicularUnit();
    *z = cross((Vector3)x, *y);
}


class Matrix4 {
public:
    Float m[4][4];
    inline Matrix4 () {memset(m, 0, sizeof m);}
    inline Matrix4 (Float diag) {
        memset(m, 0, sizeof m);
        m[0][0] = m[1][1] = m[2][2] = m[3][3] = diag;
    }
    explicit Matrix4 (const Float**from) ;
    Matrix4 (const Matrix4 & b) ;
    Matrix4 & operator = (const Matrix4 & b) ;
    const Float* operator [] (int i) const ;
    Float* operator [] (int i) ;
    Matrix4 operator * (const Matrix4 & b) const ;
};

bool inverse (const Matrix4 & _m, Matrix4 & mr) ;

class Ray {
public:
    Point3 o;
    Vector3 d;
    mutable Float t;
    inline Ray(): o(0, 0, 0), d(1, 0, 0), t(Infinity) {}
    inline Ray(const Point3 & o, const Vector3 & d, const Float & t = Infinity):
        o(o), d(d), t(t) {
    }
    inline Point3 operator () (Float dt) const {
        return o+d*dt;
    }
    inline static Ray fromTo (const Point3 & s, const Point3 & t) {
        return Ray(s, t-s, 1);
    }
};
class RayDifferential: public Ray {
public:
    bool differentials;
    Vector3 dxdir, dydir;
    Point3 dxori, dyori;
    inline RayDifferential () {}
    inline RayDifferential (const Ray & r): Ray::Ray(r) { }
    inline Ray exact (const Point2 & u) const {
        return Ray(o+((dxori-o)*u.x+(dyori-o)*u.y), (d+(dxdir-d)*u.x+(dydir-d)*u.y).normalized(), t);
    }
} ;


class Bound2 {
public:
    Point2 mn, mx;
    Bound2 ();
    explicit Bound2 (const Point2 & v);
    Bound2 (const Point2 & a, const Point2 & b) ;
};
class Bound3 {
public:
    Point3 mn, mx;
    Bound3 ();
    explicit Bound3 (const Point3 & v);
    Bound3 (const Point3 & a, const Point3 & b) ;
    bool bound (const Point3 & p) const ;
    inline Float interact(const Ray & r) const {
        if(r.o.x < mn.x && r.d.x < 0) return Infinity;
        if(r.o.x > mx.x && r.d.x > 0) return Infinity;
        if(r.o.y < mn.y && r.d.y < 0) return Infinity;
        if(r.o.y > mx.y && r.d.y > 0) return Infinity;
        if(r.o.z < mn.z && r.d.z < 0) return Infinity;
        if(r.o.z > mx.z && r.d.z > 0) return Infinity;
        if(bound(r.o)) return 0;
        Float ret = Infinity;
        if(r.d.x > 0) {
            if(mn.x >= r.o.x) {
                EFloat t = (EFloat(mn.x)-EFloat(r.o.x))/EFloat(r.d.x);
                EFloat ex = EFloat(r.o.x)+t*EFloat(r.d.x);
                EFloat ey = EFloat(r.o.y)+t*EFloat(r.d.y);
                EFloat ez = EFloat(r.o.z)+t*EFloat(r.d.z);
                if(!(ex.lower() > mx.x || ex.upper() < mn.x
                || ey.lower() > mx.y || ey.upper() < mn.y
                || ez.lower() > mx.z || ez.upper() < mn.z))
                    ret = min(ret, t.value);
            }
        } else if(r.d.x != 0) {
            if(mx.x <= r.o.x) {
                EFloat t = (EFloat(mx.x)-(EFloat)r.o.x)/EFloat(r.d.x);
                EFloat ex = EFloat(r.o.x)+t*EFloat(r.d.x);
                EFloat ey = EFloat(r.o.y)+t*EFloat(r.d.y);
                EFloat ez = EFloat(r.o.z)+t*EFloat(r.d.z);
                if(!(ex.lower() > mx.x || ex.upper() < mn.x
                || ey.lower() > mx.y || ey.upper() < mn.y
                || ez.lower() > mx.z || ez.upper() < mn.z))
                    ret = min(ret, t.value);
            }
        }
        if(r.d.y > 0) {
            if(mn.y >= r.o.y) {
                EFloat t = (EFloat(mn.y)-EFloat(r.o.y))/EFloat(r.d.y);
                EFloat ex = EFloat(r.o.x)+t*EFloat(r.d.x);
                EFloat ey = EFloat(r.o.y)+t*EFloat(r.d.y);
                EFloat ez = EFloat(r.o.z)+t*EFloat(r.d.z);
                if(!(ex.lower() > mx.x || ex.upper() < mn.x
                || ey.lower() > mx.y || ey.upper() < mn.y
                || ez.lower() > mx.z || ez.upper() < mn.z))
                    ret = min(ret, t.value);
            }
        } else if(r.d.y != 0) {
            if(mx.y <= r.o.y) {
                EFloat t = (EFloat(mx.y)-EFloat(r.o.y))/EFloat(r.d.y);
                EFloat ex = EFloat(r.o.x)+t*EFloat(r.d.x);
                EFloat ey = EFloat(r.o.y)+t*EFloat(r.d.y);
                EFloat ez = EFloat(r.o.z)+t*EFloat(r.d.z);
                if(!(ex.lower() > mx.x || ex.upper() < mn.x
                || ey.lower() > mx.y || ey.upper() < mn.y
                || ez.lower() > mx.z || ez.upper() < mn.z))
                    ret = min(ret, t.value);
            }
        }
        if(r.d.z > 0) {
            if(mn.z >= r.o.z) {
                EFloat t = (EFloat(mn.z)-EFloat(r.o.z))/EFloat(r.d.z);
                EFloat ex = EFloat(r.o.x)+t*EFloat(r.d.x);
                EFloat ey = EFloat(r.o.y)+t*EFloat(r.d.y);
                EFloat ez = EFloat(r.o.z)+t*EFloat(r.d.z);
                if(!(ex.lower() > mx.x || ex.upper() < mn.x
                || ey.lower() > mx.y || ey.upper() < mn.y
                || ez.lower() > mx.z || ez.upper() < mn.z))
                    ret = min(ret, t.value);
            }
        } else if(r.d.z != 0) {
            if(mx.z <= r.o.z) {
                EFloat t = (EFloat(mx.z)-EFloat(r.o.z))/EFloat(r.d.z);
                EFloat ex = EFloat(r.o.x)+t*EFloat(r.d.x);
                EFloat ey = EFloat(r.o.y)+t*EFloat(r.d.y);
                EFloat ez = EFloat(r.o.z)+t*EFloat(r.d.z);
                if(!(ex.lower() > mx.x || ex.upper() < mn.x
                || ey.lower() > mx.y || ey.upper() < mn.y
                || ez.lower() > mx.z || ez.upper() < mn.z))
                    ret = min(ret, t.value);
            }
        }
        return ret;
    }

    inline Float minEulerDistSquared (const Point3 & p) {
        Float dx, dy, dz;
        if(mn.x <= p.x && mx.x >= p.x) dx = 0;
        else dx = min(abs(mn.x-p.x), abs(mx.x-p.x));
        if(mn.y <= p.y && mx.y >= p.y) dy = 0;
        else dy = min(abs(mn.y-p.y), abs(mx.y-p.y));
        if(mn.z <= p.z && mx.z >= p.z) dz = 0;
        else dz = min(abs(mn.z-p.z), abs(mx.z-p.z));
        return dx*dx+dy*dy+dz*dz;
    }
};

Bound2 combine (const Bound2 & a, const Bound2 & b) ;
Bound3 combine (const Bound3 & a, const Bound3 & b) ;


class Transform {
public:
    Matrix4 m, minv;
    inline Transform () {}
    inline Transform (const Matrix4 & mt) {
        m = mt;
        inverse(m, minv);
    }
    inline Transform (const Matrix4 & mt, const Matrix4 & mtinv) {
        m = mt;
        minv = mtinv;
    }

    inline Transform operator * (const Transform & tr) const  {
        return Transform(m*tr.m, tr.minv*minv);
    }

    // Return the scales along X/Y/Z axises.
    inline Float* scales () const {
        throw "function scales todo";
    }

    inline Point3 operator () (const Point3 & v) const {
        Float t = v.x*m[3][0]+v.y*m[3][1]+v.z*m[3][2]+m[3][3];
        return Point3 (
            (v.x*m[0][0]+v.y*m[0][1]+v.z*m[0][2]+m[0][3])/t,
            (v.x*m[1][0]+v.y*m[1][1]+v.z*m[1][2]+m[1][3])/t,
            (v.x*m[2][0]+v.y*m[2][1]+v.z*m[2][2]+m[2][3])/t
        );
    }
    inline Point3 operator () (const Point3 & v, Vector3 * err) const {
        err->x = gamma(3)*(abs(v.x*m[0][0])+abs(v.y*m[0][1])+abs(v.z*m[0][2])+abs(m[0][3]));
        err->y = gamma(3)*(abs(v.x*m[1][0])+abs(v.y*m[1][1])+abs(v.z*m[1][2])+abs(m[1][3]));
        err->z = gamma(3)*(abs(v.x*m[2][0])+abs(v.y*m[2][1])+abs(v.z*m[2][2])+abs(m[2][3]));
        Float t = v.x*m[3][0]+v.y*m[3][1]+v.z*m[3][2]+m[3][3];
        return Point3 (
            ((v.x*m[0][0]+v.y*m[0][1])+(v.z*m[0][2]+m[0][3]))/t,
            ((v.x*m[1][0]+v.y*m[1][1])+(v.z*m[1][2]+m[1][3]))/t,
            ((v.x*m[2][0]+v.y*m[2][1])+(v.z*m[2][2]+m[2][3]))/t
        ) ;
    }
    inline Point3 operator () (const Point3 & v, const Vector3 & verr, Vector3 * err) const {
        err->x = 
            gamma(3)*(abs(v.x*m[0][0])+abs(v.y*m[0][1])+abs(v.z*m[0][2])+abs(m[0][3]))
        + ((Float)1+gamma(3))*(abs(verr.x*m[0][0])+abs(verr.y*m[0][1])+abs(verr.z*m[0][2]));
        err->y =
            gamma(3)*(abs(v.x*m[1][0])+abs(v.y*m[1][1])+abs(v.z*m[1][2])+abs(m[1][3]))
        + ((Float)1+gamma(3))*(abs(verr.x*m[1][0])+abs(verr.y*m[1][1])+abs(verr.z*m[1][2]));
        err->z = 
            gamma(3)*(abs(v.x*m[2][0])+abs(v.y*m[2][1])+abs(v.z*m[2][2])+abs(m[2][3]))
        + ((Float)1+gamma(3))*(abs(verr.x*m[2][0])+abs(verr.y*m[2][1])+abs(verr.z*m[2][2]));
        Float t = v.x*m[3][0]+v.y*m[3][1]+v.z*m[3][2]+m[3][3];
        return Point3 (
            ((v.x*m[0][0]+v.y*m[0][1])+(v.z*m[0][2]+m[0][3]))/t,
            ((v.x*m[1][0]+v.y*m[1][1])+(v.z*m[1][2]+m[1][3]))/t,
            ((v.x*m[2][0]+v.y*m[2][1])+(v.z*m[2][2]+m[2][3]))/t
        ) ;
    }
    inline Vector3 operator () (const Vector3 & v) const {
        return Vector3(
            v.x*m[0][0]+v.y*m[0][1]+v.z*m[0][2],
            v.x*m[1][0]+v.y*m[1][1]+v.z*m[1][2],
            v.x*m[2][0]+v.y*m[2][1]+v.z*m[2][2]
        );
    }
    inline Vector3 operator () (const Vector3 & v, Vector3 * err) const {
        err->x = gamma(3)*(abs(v.x*m[0][0])+abs(v.y*m[0][1])+abs(v.z*m[0][2])+abs(m[0][3]));
        err->y = gamma(3)*(abs(v.x*m[1][0])+abs(v.y*m[1][1])+abs(v.z*m[1][2])+abs(m[1][3]));
        err->z = gamma(3)*(abs(v.x*m[2][0])+abs(v.y*m[2][1])+abs(v.z*m[2][2])+abs(m[2][3]));
        return Vector3(
            v.x*m[0][0]+v.y*m[0][1]+v.z*m[0][2],
            v.x*m[1][0]+v.y*m[1][1]+v.z*m[1][2],
            v.x*m[2][0]+v.y*m[2][1]+v.z*m[2][2]
        );
    }
    inline Normal operator () (const Normal & v) const {
        return Normal(Vector3 (
            v.x*minv[0][0]+v.y*minv[1][0]+v.z*minv[2][0],
            v.x*minv[0][1]+v.y*minv[1][1]+v.z*minv[2][1],
            v.x*minv[0][2]+v.y*minv[1][2]+v.z*minv[2][2]
        ).normalizedSafe());
    }

    // will offset the ray origin automatically
    inline Ray operator () (const Ray &v) const {
        Vector3 err;
        Point3 no = (*this)(v.o, &err);
        Vector3 nd = (*this)(v.d);
        Float ndl = dot(abs(nd), err)/std::sqrt(nd.lengthSquared());
        return Ray(no+nd*ndl, nd, v.t-ndl);
    }
    // will NOT offset the ray origin automatically
    inline Ray operator () (const Ray &v, Vector3 * erro) const {
        Point3 no = (*this)(v.o, erro);
        Vector3 nd = (*this)(v.d);
        return Ray(no, nd, v.t);
    }
    // will NOT offset the ray origin automatically
    inline Ray operator () (const Ray &v, Vector3 * erro, Vector3 * errd) const {
        Point3 no = (*this)(v.o, erro);
        Vector3 nd = (*this)(v.d, errd);
        return Ray(no, nd, v.t);
    }

    inline Bound3 operator () (const Bound3 & bnd) const 
    {
        Point3 p1 = (*this)(Point3(bnd.mn.x, bnd.mn.y, bnd.mn.z));
        Point3 p2 = (*this)(Point3(bnd.mn.x, bnd.mn.y, bnd.mx.z));
        Point3 p3 = (*this)(Point3(bnd.mn.x, bnd.mx.y, bnd.mn.z));
        Point3 p4 = (*this)(Point3(bnd.mn.x, bnd.mx.y, bnd.mx.z));
        Point3 p5 = (*this)(Point3(bnd.mx.x, bnd.mn.y, bnd.mn.z));
        Point3 p6 = (*this)(Point3(bnd.mx.x, bnd.mn.y, bnd.mx.z));
        Point3 p7 = (*this)(Point3(bnd.mx.x, bnd.mx.y, bnd.mn.z));
        Point3 p8 = (*this)(Point3(bnd.mx.x, bnd.mx.y, bnd.mx.z));
        Float xmin = min(
            min(min(p1.x, p2.x), min(p3.x, p4.x)),
            min(min(p5.x, p6.x), min(p7.x, p8.x))
        );
        Float ymin = min(
            min(min(p1.y, p2.y), min(p3.y, p4.y)),
            min(min(p5.y, p6.y), min(p7.y, p8.y))
        );
        Float zmin = min(
            min(min(p1.z, p2.z), min(p3.z, p4.z)),
            min(min(p5.z, p6.z), min(p7.z, p8.z))
        );
        Float xmax = max(
            max(max(p1.x, p2.x), max(p3.x, p4.x)),
            max(max(p5.x, p6.x), max(p7.x, p8.x))
        );
        Float ymax = max(
            max(max(p1.y, p2.y), max(p3.y, p4.y)),
            max(max(p5.y, p6.y), max(p7.y, p8.y))
        );
        Float zmax = max(
            max(max(p1.z, p2.z), max(p3.z, p4.z)),
            max(max(p5.z, p6.z), max(p7.z, p8.z))
        );
        return Bound3(Point3(xmin, ymin, zmin), Point3(xmax, ymax, zmax));
    }

    inline RayDifferential operator () (const RayDifferential & v) const 
    {
        RayDifferential nd = (*this)(v);
        nd.dxori = (*this)(v.dxori);
        nd.dyori = (*this)(v.dyori);
        nd.dxdir = (*this)(v.dxdir);
        nd.dydir = (*this)(v.dydir);
        return nd;
    }

    inline Point3 reversed (const Point3 & v) const {
        Float t = v.x*minv[3][0]+v.y*minv[3][1]+v.z*minv[3][2]+minv[3][3];
        return Point3 (
            (v.x*minv[0][0]+v.y*minv[0][1]+v.z*minv[0][2]+minv[0][3])/t,
            (v.x*minv[1][0]+v.y*minv[1][1]+v.z*minv[1][2]+minv[1][3])/t,
            (v.x*minv[2][0]+v.y*minv[2][1]+v.z*minv[2][2]+minv[2][3])/t
        );
    }
    inline Point3 reversed (const Point3 & v, Vector3 * err) const {
        err->x = gamma(3)*(abs(v.x*minv[0][0])+abs(v.y*minv[0][1])+abs(v.z*minv[0][2])+abs(minv[0][3]));
        err->y = gamma(3)*(abs(v.x*minv[1][0])+abs(v.y*minv[1][1])+abs(v.z*minv[1][2])+abs(minv[1][3]));
        err->z = gamma(3)*(abs(v.x*minv[2][0])+abs(v.y*minv[2][1])+abs(v.z*minv[2][2])+abs(minv[2][3]));
        Float t = v.x*minv[3][0]+v.y*minv[3][1]+v.z*minv[3][2]+minv[3][3];
        return Point3(
            ((v.x*minv[0][0]+v.y*minv[0][1])+(v.z*minv[0][2]+minv[0][3]))/t,
            ((v.x*minv[1][0]+v.y*minv[1][1])+(v.z*minv[1][2]+minv[1][3]))/t,
            ((v.x*minv[2][0]+v.y*minv[2][1])+(v.z*minv[2][2]+minv[2][3]))/t
        ) ;
    }
    inline Point3 reversed (const Point3 & v, const Vector3 & verr, Vector3 * err) const {
        err->x =
            gamma(3)*(abs(v.x*minv[0][0])+abs(v.y*minv[0][1])+abs(v.z*minv[0][2])+abs(minv[0][3]))
        + ((Float)1+gamma(3))*(abs(verr.x*minv[0][0])+abs(verr.y*minv[0][1])+abs(verr.z*minv[0][2]));
        err->y =
            gamma(3)*(abs(v.x*minv[1][0])+abs(v.y*minv[1][1])+abs(v.z*minv[1][2])+abs(minv[1][3]))
        + ((Float)1+gamma(3))*(abs(verr.x*minv[1][0])+abs(verr.y*minv[1][1])+abs(verr.z*minv[1][2]));
        err->z =
            gamma(3)*(abs(v.x*minv[2][0])+abs(v.y*minv[2][1])+abs(v.z*minv[2][2])+abs(minv[2][3]))
        + ((Float)1+gamma(3))*(abs(verr.x*minv[2][0])+abs(verr.y*minv[2][1])+abs(verr.z*minv[2][2]));
        Float t = v.x*minv[3][0]+v.y*minv[3][1]+v.z*minv[3][2]+minv[3][3];
        return Point3( // adjust computation order to minimize upper boundary of error
            ((v.x*minv[0][0]+v.y*minv[0][1])+(v.z*minv[0][2]+minv[0][3]))/t, // Question: why not taking t into account?
            ((v.x*minv[1][0]+v.y*minv[1][1])+(v.z*minv[1][2]+minv[1][3]))/t,
            ((v.x*minv[2][0]+v.y*minv[2][1])+(v.z*minv[2][2]+minv[2][3]))/t
        ) ; 
    }
    inline Vector3 reversed (const Vector3 & v) const {
        return Vector3(
            v.x*minv[0][0]+v.y*minv[0][1]+v.z*minv[0][2],
            v.x*minv[1][0]+v.y*minv[1][1]+v.z*minv[1][2],
            v.x*minv[2][0]+v.y*minv[2][1]+v.z*minv[2][2]
        );
    }
    inline Vector3 reversed (const Vector3 & v, Vector3 * err) const {
        err->x = gamma(4)*(abs(v.x*minv[0][0])+abs(v.y*minv[0][1])+abs(v.z*minv[0][2]));
        err->y = gamma(4)*(abs(v.x*minv[1][0])+abs(v.y*minv[1][1])+abs(v.z*minv[1][2]));
        err->z = gamma(4)*(abs(v.x*minv[2][0])+abs(v.y*minv[2][1])+abs(v.z*minv[2][2]));
        return Vector3(
            v.x*minv[0][0]+v.y*minv[0][1]+v.z*minv[0][2],
            v.x*minv[1][0]+v.y*minv[1][1]+v.z*minv[1][2],
            v.x*minv[2][0]+v.y*minv[2][1]+v.z*minv[2][2]
        );
    }
    inline Normal reversed (const Normal & v) const {
        return Normal(Vector3(
            v.x*m[0][0]+v.y*m[1][0]+v.z*m[2][0],
            v.x*m[0][1]+v.y*m[1][1]+v.z*m[2][1],
            v.x*m[0][2]+v.y*m[1][2]+v.z*m[2][2]
        ).normalizedSafe());
    }

    // will offset the ray origin automatically
    inline Ray reversed (const Ray &v) const {
        Vector3 err;
        Point3 no = this->reversed(v.o, &err);
        Vector3 nd = this->reversed(v.d);
        Float ndl = dot(abs(nd), err)/std::sqrt(nd.lengthSquared());
        return Ray((Point3)(no+nd*ndl), nd, v.t-ndl);
    }

    // will NOT offset the ray origin automatically
    Ray reversed (const Ray &v, Vector3 * erro) const { 
        Point3 no = this->reversed(v.o, erro);
        Vector3 nd = this->reversed(v.d);
        return Ray(no, nd, v.t);
    }

    // will NOT offset the ray origin automatically
    inline Ray reversed (const Ray &v, Vector3 * erro, Vector3 * errd) const { 
        Point3 no = this->reversed(v.o, erro);
        Vector3 nd = this->reversed(v.d, errd);
        return Ray(no, nd, v.t);
    }

    inline Bound3 reversed (const Bound3 & bnd) const 
    {

        Point3 p1 = reversed(Point3(bnd.mn.x, bnd.mn.y, bnd.mn.z));
        Point3 p3 = reversed(Point3(bnd.mn.x, bnd.mx.y, bnd.mn.z));
        Point3 p2 = reversed(Point3(bnd.mn.x, bnd.mn.y, bnd.mx.z));
        Point3 p4 = reversed(Point3(bnd.mn.x, bnd.mx.y, bnd.mx.z));
        Point3 p5 = reversed(Point3(bnd.mx.x, bnd.mn.y, bnd.mn.z));
        Point3 p6 = reversed(Point3(bnd.mx.x, bnd.mn.y, bnd.mx.z));
        Point3 p7 = reversed(Point3(bnd.mx.x, bnd.mx.y, bnd.mn.z));
        Point3 p8 = reversed(Point3(bnd.mx.x, bnd.mx.y, bnd.mx.z));
        Float xmin = min(
            min(min(p1.x, p2.x), min(p3.x, p4.x)),
            min(min(p5.x, p6.x), min(p7.x, p8.x))
        );
        Float ymin = min(
            min(min(p1.y, p2.y), min(p3.y, p4.y)),
            min(min(p5.y, p6.y), min(p7.y, p8.y))
        );
        Float zmin = min(
            min(min(p1.z, p2.z), min(p3.z, p4.z)),
            min(min(p5.z, p6.z), min(p7.z, p8.z))
        );
        Float xmax = max(
            max(max(p1.x, p2.x), max(p3.x, p4.x)),
            max(max(p5.x, p6.x), max(p7.x, p8.x))
        );
        Float ymax = max(
            max(max(p1.y, p2.y), max(p3.y, p4.y)),
            max(max(p5.y, p6.y), max(p7.y, p8.y))
        );
        Float zmax = max(
            max(max(p1.z, p2.z), max(p3.z, p4.z)),
            max(max(p5.z, p6.z), max(p7.z, p8.z))
        );
        return Bound3(Point3(xmin, ymin, zmin), Point3(xmax, ymax, zmax));    }

    inline Transform reversed() const 
    {
        return Transform(minv, m);
    }

    inline RayDifferential reversed(const RayDifferential & v) const 
    {
        RayDifferential nd;
        nd.dxori = reversed(v.dxori);
        nd.dyori = reversed(v.dyori);
        nd.dxdir = reversed(v.dxdir);
        nd.dydir = reversed(v.dydir);
        return nd;
    }

    Interact operator () (const Interact & in) const ;
    Interact reversed (const Interact & in) const ;
    SurfInteract operator () (const SurfInteract & si) const ;
    SurfInteract reversed (const SurfInteract & si) const ;
};
extern Transform transform_unit;
Transform scale (Vector3 s) ;
Transform shift (Vector3 v) ;
Transform rotateX (Float f) ;
Transform rotateY (Float f) ;
Transform rotateZ (Float f) ;


}
#endif
